#Loading packages
library(jtools)
library(sandwich)
library(sjPlot)
library(sjlabelled)
library(sjmisc)
library(ggplot2)
library(car)
library(lmtest)
library(multiwayvcov)
library(plyr)
library(dummies)
library(gridExtra)
library(dplyr)
library(tidyr)
library(stargazer)
library(broom)
library(readxl)
library(tidyverse)
library(readstata13)
library(haven)
library(estimatr)
library(magrittr) 
library(texreg)
library(xtable)
library(gridExtra)
library(ggExtra)
library(cobalt)
library(MatchIt)
library(WeightIt)
library(Rmisc)
library(reshape2)
library(base)
library(dotwhisker)
library(broom)
library(tidyverse)
library(hrbrthemes)
library(viridis)
library(forcats)
library(scales)
library(readxl)
library(xlsx)
library(openxlsx)
library(readstata13)
library(devtools)
library(factoextra)
library(cregg)
library(dplyr)

## Please set your directory here
councilors <- read.xlsx("councilors_profiledata_withDV2_PCA_weight.xlsx")

councilors$Proximity <- factor(councilors$Proximity)
councilors$Public_goods <- factor(councilors$Publicgoods)
councilors$Size <- factor(councilors$Size)
councilors$Type <- factor(councilors$Type)
councilors$Run_by <- factor(councilors$Runby)
councilors$munic_res <- factor(councilors$munic_res)
dv_binary <- factor(councilors$dv_binary)
councilors$ResponseId <- c(councilors$ResponseId)

attr(councilors$Run_by, "label") <- "Run by"
attr(councilors$Public_goods, "label") <- "Public goods"

councilors$Size <- factor(councilors$Size, levels =
                            c("< 1% of local population","1% of local population", "> 1% of local population"))
councilors$Proximity <- factor(councilors$Proximity, levels =
                                 c("In ctr", "< 30mins from ctr", "> 30mins from ctr"))
councilors$Type <- factor(councilors$Type, levels =
                            c("Fully open", "Partially open", "Closed"))

# Figure 1:  Aggregate marginal means (MMs)

f1 <- dv_binary ~ Proximity + Public_goods + Size+ Type + Run_by
plot(mm(councilors, f1, id = ~ResponseId), vline = 0.5)

d1 <- mm(councilors, dv_binary ~ Proximity + Public_goods + Size+ Type + Run_by,
         id = ~ ResponseId)


apatheme=theme_bw()+
  theme(panel.grid.major=element_blank(),
        axis.line=element_line(),
        text=element_text(family='Helvetica'),
        axis.text=element_text(size=11),
        axis.title=element_text(size=11),
        legend.text = element_text(size = 11),
        legend.title = element_text(size=14))

apatheme2=theme_bw()+
  theme(panel.grid.major=element_blank(),
        axis.line=element_line(),
        text=element_text(family='Helvetica'),
        axis.text=element_text(size=11),
        axis.title=element_text(size=11),
        legend.text = element_text(size = 11),
        legend.title = element_blank())

F1 <- plot(d1, feature_headers = FALSE, vline = 0.5) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right")

png("F1.png", width = 600, height = 650)
F1 + apatheme2
dev.off()


# Figure 2:  MMs based on presence of active host sites
f3 <- dv_binary ~ Size
councilors$treat1 <- as.numeric(councilors$treat1)
table(councilors$treat1)
councilors$Camp[councilors$treat1==1]<-"Yes"
councilors$Camp[councilors$treat1==0]<-"No"

councilors$Camp <- factor(councilors$Camp,levels = c("Yes", "No"))

cj(councilors, f3, id = ~ResponseId, by = ~Camp, estimate = "mm")

F2 <- plot(cj(councilors, f3, id = ~ResponseId, by = ~Camp, estimate = "mm"),
           group = "Camp", vline = 0.5)

F2 <- F2+apatheme
png("F2.png", width = 400, height = 400)
F2
dev.off()

# Figure B.2:  Aggregate marginal means showing profile placement diagnostics
councilors$profile_number <- factor(councilors$profile_number)

councilors$Profile[councilors$profile_number==1]<-"Left-profile"
councilors$Profile[councilors$profile_number==2]<-"Right-profile"
councilors$Profile <- factor(councilors$Profile, levels= c("Left-profile", "Right-profile"))
mm <- cj(councilors, f1, id = ~ResponseId, estimate = "mm", 
         by = ~Profile)
B2 <- plot(rbind(mm),vline = 0.5) + ggplot2::facet_wrap(~BY, ncol = 2L)
B2 <- B2 + apatheme

B2

cj_freqs(councilors, f1, id = ~ResponseId)

# Figure B.3:  Aggregate marginal means with Likert scale outcome

councilors$dv_scale22 <- (councilors$dv_scale-0)/7
summary(councilors$dv_scale22)

f4 <- dv_scale22 ~ Proximity + Public_goods + Size+ Type + Run_by
B3 <- plot(mm(councilors, f4, id = ~ResponseId))

png("B3.png", width = 600, height = 650)
B3+apatheme
dev.off()

# Figure B.4:  Aggregate marginal means with entropy balancing weights
webal <- councilors$webal
webal2 <- councilors$webal2

w1 <- plot(mm(councilors, f1,
              id = ~ResponseId,
              weights = ~webal), vline = 0.5)
w1 <- w1 + ggtitle("Weight 1 (gender, geog.)")

w2 <- plot(mm(councilors, f1,
              id = ~ResponseId,
              weights = ~webal2), vline = 0.5)
w2 <- w2 + ggtitle("Weight 2 (gender, geog., existing camp)")

grid.arrange(w1+ apatheme, w2 + apatheme, nrow=1)


# Table C.2:  Aggregate average marginal component effects (AMCEs)
library(fastDummies)

councilors <- dummy_cols(councilors, select_columns = "Proximity")
councilors <- dummy_cols(councilors, select_columns = "Public_goods")
councilors <- dummy_cols(councilors, select_columns = "Size")
councilors <- dummy_cols(councilors, select_columns = "Type")
councilors <- dummy_cols(councilors, select_columns = "Run_by")


lessthan1pct <- c(councilors$`Size_< 1% of local population`)
onepct <- c(councilors$`Size_1% of local population`)
morethan1pct <- c(councilors$`Size_> 1% of local population`)
lessthan30minwalk <- c(councilors$`Proximity_< 30mins from ctr`)
morethan30minwalk <- c(councilors$`Proximity_> 30mins from ctr`)
inthecenter <- c(councilors$`Proximity_In ctr`)
army <- c(councilors$Run_by_Army)
church <- c(councilors$Run_by_Church)
government <- c(councilors$Run_by_Government)
IOs <- c(councilors$`Run_by_IOs (UNHCR,IOM)`)
local_government <- c(councilors$`Run_by_Local Government`)
hire_more_municipal_employees <- c(councilors$`Public_goods_Hire more municipal employees`)
hire_more_teachers_and_doctors <- c(councilors$`Public_goods_Hire more teachers and doctors`)
more_infrastructure <- c(councilors$`Public_goods_More infrastructure to the municipality`)
closed <- c(councilors$Type_Closed)
partially_open <- c(councilors$`Type_Partially open`)
fully_open <- c(councilors$`Type_Fully open`)
ResponseId <- c(councilors$ResponseId)

reg1 <- lm_robust(dv_binary ~ morethan30minwalk + lessthan30minwalk
                  + more_infrastructure + hire_more_teachers_and_doctors
                  + morethan1pct +morethan1pct+ lessthan1pct
                  + partially_open + fully_open +
                    church + government + IOs + local_government,
                  clusters = ResponseId, data=councilors)

reg2 <- lm_robust(dv_binary ~ morethan30minwalk + lessthan30minwalk
                  + more_infrastructure + hire_more_teachers_and_doctors
                  + morethan1pct +morethan1pct+ lessthan1pct
                  + partially_open + fully_open +
                    church + government + IOs + local_government,
                  clusters = ResponseId, data=councilors, weights = webal2)

reg3 <- lm_robust(dv_binary ~ morethan30minwalk + lessthan30minwalk
                  + more_infrastructure + hire_more_teachers_and_doctors
                  + morethan1pct +morethan1pct+ lessthan1pct
                  + partially_open + fully_open +
                    church + government + IOs + local_government,
                  clusters = ResponseId, data=councilors,
                  fixed_effects = munic_res)

texreg(list(reg1, reg2, reg3),
       custom.model.names = c("Main model", "Weighted model",
                              "Municipality FE model"), include.ci=FALSE)


# Table C.3:  Aggregate marginal means (MMs)
agg <- mm(councilors, f1,
   id = ~ResponseId)
print(xtable(agg, type = "latex"), file = "mm_aggregate.tex")


# Table C.4:  Active host site marginal means (MMs)
ac <- cj(councilors, f3, id = ~ResponseId, by = ~Camp, estimate = "mm")
print(xtable(ac, type = "latex"), file = "mm_hostsite.tex")

# Table C.5: Subgroup analysis, differences in marginal means
test<-mm_diffs(councilors, f3, by = ~ Camp, id = ~ ResponseId)
print(xtable(test, type = "latex"), file = "mm_diff_test.tex")

# H_null: mm_camp=mm_nocamp, H_alt: mm_camp<mm_nocamp => mm_nocamp-mm_camp>0
# critical value z from standard normal distribution is 1.64 for ??=5%.
# -0.9877<critical value z, therefore we fail to reject the null hypothesis 

#create p-value, P(Z>z)=P(Z>-0.9877)=1-P(Z=<-0.9877)=1-P(Z>0.9877)=P(Z=<0.9877)
pvalue<-pnorm(0.9877, lower.tail = TRUE)

# Table C.6:  Active host site AMCEs
treated_group = subset(councilors, councilors$treat1==1)

lessthan1pct <- c(treated_group$`Size_< 1% of local population`)
onepct <- c(treated_group$`Size_1% of local population`)
morethan1pct <- c(treated_group$`Size_> 1% of local population`)
lessthan30minwalk <- c(treated_group$`Proximity_< 30mins from ctr`)
morethan30minwalk <- c(treated_group$`Proximity_> 30mins from ctr`)
inthecenter <- c(treated_group$`Proximity_In ctr`)
army <- c(treated_group$Run_by_Army)
church <- c(treated_group$Run_by_Church)
government <- c(treated_group$Run_by_Government)
IOs <- c(treated_group$`Run_by_IOs (UNHCR,IOM)`)
local_government <- c(treated_group$`Run_by_Local Government`)
hire_more_municipal_employees <- c(treated_group$`Public_goods_Hire more municipal employees`)
hire_more_teachers_and_doctors <- c(treated_group$`Public_goods_Hire more teachers and doctors`)
more_infrastructure <- c(treated_group$`Public_goods_More infrastructure to the municipality`)
closed <- c(treated_group$Type_Closed)
partially_open <- c(treated_group$`Type_Partially open`)
fully_open <- c(treated_group$`Type_Fully open`)
ResponseId <- c(treated_group$ResponseId)

reg1_councilors1 <- lm_robust(dv_binary ~ morethan30minwalk + lessthan30minwalk
                              + more_infrastructure + hire_more_teachers_and_doctors
                              + morethan1pct +morethan1pct+ lessthan1pct
                              + partially_open + fully_open +
                                church + government + IOs + local_government,
                              clusters = ResponseId, data=treated_group)
summary(reg1_councilors1)

nontreated_group = subset(councilors, councilors$treat1==0)

lessthan1pct <- c(nontreated_group$`Size_< 1% of local population`)
onepct <- c(nontreated_group$`Size_1% of local population`)
morethan1pct <- c(nontreated_group$`Size_> 1% of local population`)
lessthan30minwalk <- c(nontreated_group$`Proximity_< 30mins from ctr`)
morethan30minwalk <- c(nontreated_group$`Proximity_> 30mins from ctr`)
inthecenter <- c(nontreated_group$`Proximity_In ctr`)
army <- c(nontreated_group$Run_by_Army)
church <- c(nontreated_group$Run_by_Church)
government <- c(nontreated_group$Run_by_Government)
IOs <- c(nontreated_group$`Run_by_IOs (UNHCR,IOM)`)
local_government <- c(nontreated_group$`Run_by_Local Government`)
hire_more_municipal_employees <- c(nontreated_group$`Public_goods_Hire more municipal employees`)
hire_more_teachers_and_doctors <- c(nontreated_group$`Public_goods_Hire more teachers and doctors`)
more_infrastructure <- c(nontreated_group$`Public_goods_More infrastructure to the municipality`)
closed <- c(nontreated_group$Type_Closed)
partially_open <- c(nontreated_group$`Type_Partially open`)
fully_open <- c(nontreated_group$`Type_Fully open`)
ResponseId <- c(nontreated_group$ResponseId)

reg1_councilors2 <- lm_robust(dv_binary ~ morethan30minwalk + lessthan30minwalk
                              + more_infrastructure + hire_more_teachers_and_doctors
                              + morethan1pct +morethan1pct+ lessthan1pct
                              + partially_open + fully_open +
                                church + government + IOs + local_government,
                              clusters = ResponseId, data=nontreated_group)
summary(reg1_councilors2)

texreg(list(reg1_councilors1, reg1_councilors2),
       custom.model.names = c("Treated", "Not treated"),
       include.ci = FALSE)

# Table D.7:  Descriptive statistics
summary(councilors$cultural_threat)
summary(councilors$economic_threat)
mean(councilors$ethnocentric_pca, na.rm=TRUE)
councilors$ethnocentric_pca[is.na(councilors$ethnocentric_pca)]<-6.415777
summary(councilors$ethnocentric_pca)
summary(councilors$Q16_1)

stargazer(councilors[c("cultural_threat","economic_threat","ethnocentric_pca", "Q16_1")],
          title="Descriptive statistics", digits=3)

# Figure D.6:  Histogram of self-reported ideology

ideology <- as.numeric(councilors$Q16_1)
hist(ideology,  col="cadetblue")

# Table D.12:  Subgroup AMCEs

# PCA sociocultural threat
summary(councilors$cultural_threat)
councilors_lowthreat <- subset(councilors, cultural_threat < 1.410)
councilors_highthreat <- subset(councilors, cultural_threat >= 1.410)
# low threat #
dv_binary <- c(councilors_lowthreat$dv_binary)
response_id <- c(councilors_lowthreat$ResponseId)
reg1_councilors_scl <- lm_robust(dv_binary ~ Proximity + Public_goods + Size+ Type + Run_by,
                                 clusters = response_id, data = councilors_lowthreat)
# high threat #
dv_binary <- c(councilors_highthreat$dv_binary)
response_id <- c(councilors_highthreat$ResponseId)
reg1_councilors_sch <- lm_robust(dv_binary ~ Proximity + Public_goods + Size+ Type + Run_by,
                                 clusters = response_id, data = councilors_highthreat)

# PCA economic threat #
summary(councilors$economic_threat)
councilors_lowthreat <- subset(councilors, economic_threat < 3.525)
councilors_highthreat <- subset(councilors, economic_threat >= 3.525)
# low threat #
dv_binary <- c(councilors_lowthreat$dv_binary)
response_id <- c(councilors_lowthreat$ResponseId)
reg1_councilors_el <- lm_robust(dv_binary ~ Proximity + Public_goods + Size
                                + Type + Run_by, clusters = response_id, data=councilors_lowthreat)
# high threat #
dv_binary <- c(councilors_highthreat$dv_binary)
response_id <- c(councilors_highthreat$ResponseId)
reg1_councilors_eh <- lm_robust(dv_binary ~ Proximity + Public_goods + Size
                                + Type + Run_by, clusters = response_id, data=councilors_highthreat)
# PCA ethnocentric values #
summary(councilors$ethnocentric_pca)
councilors_lowthreat <- subset(councilors, ethnocentric_pca < 6.544)
councilors_highthreat <- subset(councilors, ethnocentric_pca >= 6.544)
# low #
dv_binary <- c(councilors_lowthreat$dv_binary)
response_id <- c(councilors_lowthreat$ResponseId)
reg1_councilors_vl<- lm_robust(dv_binary ~ Proximity + Public_goods + Size
                               + Type + Run_by, clusters = response_id, data=councilors_lowthreat)
# high #
dv_binary <- c(councilors_highthreat$dv_binary)
response_id <- c(councilors_highthreat$ResponseId)
reg1_councilors_vh <- lm_robust(dv_binary ~ Proximity + Public_goods + Size
                                + Type + Run_by, clusters = response_id,data=councilors_highthreat)

# Ideology #
ideology <- as.numeric(councilors$Q16_1)
councilors_right <- subset(councilors, ideology>=5)
councilors_left <- subset(councilors, ideology<5)
# right #
dv_binary <- c(councilors_right$dv_binary)
response_id <- c(councilors_right$ResponseId)
reg1_councilors_right <- lm_robust(dv_binary ~ Proximity + Public_goods + Size
                                   + Type + Run_by, clusters = response_id, data=councilors_right)
# left #
dv_binary <- c(councilors_left$dv_binary)
response_id <- c(councilors_left$ResponseId)
reg1_councilors_left <- lm_robust(dv_binary ~ Proximity + Public_goods + Size
                                  + Type + Run_by, clusters = response_id, data=councilors_left)

texreg(list(reg1_councilors_scl, reg1_councilors_sch, reg1_councilors_el, reg1_councilors_eh, reg1_councilors_vl, reg1_councilors_vh, reg1_councilors_left, reg1_councilors_right),
       custom.model.names = c("SC-low", "SC-high", "E-low", "E-high", "V-low", "V-high", "Left", "Right"),
       include.ci=FALSE)

# Figure D.7: Subgroup marginal means - PCAs (MMs)

apatheme3=theme_bw()+
  theme(panel.grid.major=element_blank(),
        axis.line=element_line(),
        text=element_text(family='Helvetica'),
        axis.text=element_text(size=10),
        axis.title=element_text(size=10),
        legend.text = element_text(size = 10),
        legend.title = element_text(size=13))

# sociocultural threat 
summary(councilors$cultural_threat)
councilors <- councilors %>% mutate(Sociocultural_threat = case_when(cultural_threat >= 1.410  ~ 'High',
                                                                     cultural_threat < 1.410  ~ 'Low'))

councilors$Sociocultural_threat <- factor(councilors$Sociocultural_threat, levels = c('High', 'Low'))

f2 <- dv_binary ~ Type + Run_by
f1 <- dv_binary ~ Proximity + Public_goods + Size+ Type + Run_by

s1 <- plot(cj(councilors, f1, id = ~ResponseId, by = ~Sociocultural_threat, estimate = "mm"),
           group = "Sociocultural_threat", vline = 0.5)

s11 <- plot(s1, feature_headers = FALSE, vline = 0.5) + apatheme3

# economic threat #
summary(councilors$economic_threat)
councilors <- councilors %>% mutate(Economic_threat = case_when(economic_threat >= 3.525  ~ 'High',
                                                                economic_threat < 3.525 ~ 'Low'))
councilors$Economic_threat <- factor(councilors$Economic_threat, levels = c('High', 'Low'))

f2 <- dv_binary ~ Type + Run_by
f1 <- dv_binary ~ Proximity + Public_goods + Size+ Type + Run_by

s2 <- plot(cj(councilors, f1, id = ~ResponseId, by = ~Economic_threat, estimate = "mm"), group = "Economic_threat", vline = 0.5)

s22 <- plot(s2, feature_headers = FALSE, vline = 0.5) + apatheme3


# ethnocentric values #
summary(councilors$ethnocentric_pca)
councilors <- councilors %>% mutate(Ethnocentric_values = case_when(ethnocentric_pca >= 6.544  ~ 'High',
                                                                    ethnocentric_pca < 6.544  ~ 'Low'))
councilors$Ethnocentric_values <- factor(councilors$Ethnocentric_values, levels = c('High', "Low"))

f2 <- dv_binary ~ Type + Run_by
f1 <- dv_binary ~ Proximity + Public_goods + Size+ Type + Run_by

s3 <- plot(cj(councilors, f1, id = ~ResponseId, by = ~Ethnocentric_values, estimate = "mm"), group = "Ethnocentric_values", vline = 0.5)

s33 <- plot(s3, feature_headers = FALSE, vline = 0.5) + apatheme3
grid.arrange(s11,s22,s33, ncol = 3, layout_matrix=rbind(c(1,1,2,2), c(NA, 3, 3, NA)))

# Figure D.8: Subgroup marginal means - ideology (MMs)
ideology <- as.numeric(councilors$Q16_1)
summary(ideology)
councilors <- councilors %>% mutate(Reported_Ideology = case_when(ideology >= 5  ~ 'Right',
                                                                  ideology < 5  ~ 'Left'))

councilors$Reported_Ideology <- factor(councilors$Reported_Ideology, levels = c('Right', 'Left'))

f2 <- dv_binary ~ Type + Run_by
f1 <- dv_binary ~ Proximity + Public_goods + Size+ Type + Run_by

s4 <- plot(cj(councilors, f1, id = ~ResponseId, by = ~Reported_Ideology, estimate = "mm"),
           group = "Reported_Ideology", vline = 0.5)

s44 <- plot(s4, feature_headers = FALSE, vline = 0.5) + apatheme
s44



# Table D.13:  Subgroup marginal means (MMs)
a <- cj(councilors, f1, id = ~ResponseId, by = ~Sociocultural_threat, estimate = "mm")
b <- cj(councilors, f1, id = ~ResponseId, by = ~Economic_threat, estimate = "mm")
c <- cj(councilors, f1, id = ~ResponseId, by = ~Ethnocentric_values, estimate = "mm")
d <- cj(councilors, f1, id = ~ResponseId, by = ~Reported_Ideology, estimate = "mm")

print(xtable(a, type = "latex"), file = "mm1.tex")
print(xtable(b, type = "latex"), file = "mm2.tex")
print(xtable(c, type = "latex"), file = "mm3.tex")
print(xtable(d, type = "latex"), file = "mm4.tex")


# Figure D.9: Ideology sensitivity checks (MMs)
f2 <- dv_binary ~ Type + Run_by
f1 <- dv_binary ~ Proximity + Public_goods + Size+ Type + Run_by

# mid-point cutoff
ideology <- as.numeric(councilors$Q16_1)
summary(ideology)
councilors <- councilors %>% mutate(Ideology = case_when(ideology >= 5  ~ 'Right',
                                                         ideology < 5  ~ 'Left'))

councilors$Ideology <- factor(councilors$Ideology, levels = c('Right', 'Left'))
u <- plot(cj(councilors, f1, id = ~ResponseId, by = ~Ideology, estimate = "mm"),
          group = "Ideology", vline = 0.5)

u <- u + ggtitle("Median-point cut-off")

# median + 1 -1 cutoff
ideology <- as.numeric(councilors$Q16_1)
summary(ideology)

councilors <- councilors %>% mutate(Ideology = case_when(ideology >= 6  ~ 'Right (+1)',
                                                         ideology < 4  ~ 'Left (-1)'))

councilors$Ideology <- factor(councilors$Ideology, levels = c('Right (+1)', 'Left (-1)'))
u2 <- plot(cj(councilors, f1, id = ~ResponseId, by = ~Ideology, estimate = "mm"),
           group = "Ideology", vline = 0.5)
u2 <- u2 + ggtitle("Median-point +1/-1 cut-off")

# mean-point cutoff
ideology <- as.numeric(councilors$Q16_1)
summary(ideology)
councilors <- councilors %>% mutate(Ideology = case_when(ideology >= 5.02  ~ 'Right',
                                                         ideology < 5.02  ~ 'Left'))

councilors$Ideology <- factor(councilors$Ideology, levels = c('Right', 'Left'))
u3 <- plot(cj(councilors, f1, id = ~ResponseId, by = ~Ideology, estimate = "mm"),
           group = "Ideology", vline = 0.5)
u3 <- u3 + ggtitle("Mean-point cut-off")

# tertile-point cutoff
# Find tertiles
vTert = quantile(councilors$Q16_1, c(0:3/3), na.rm = TRUE)

# classify values
councilors$Ideology = with(councilors, 
               cut(Q16_1, 
                   vTert, 
                   include.lowest = T, 
                   labels = c("Left", "Center", "Right")))

u4 <- plot(cj(councilors, f1, id = ~ResponseId, by = ~Ideology, estimate = "mm"),
           group = "Ideology", vline = 0.5)
u4 <- u4 + ggtitle("Tertile cut-off")

grid.arrange(u+apatheme3,u2+apatheme3,u3+apatheme3,u4+apatheme3)
#save from viewer

# Figure D.10: Sociocultural threat PCA sensitivity checks
summary(councilors$cultural_threat)

councilors <- councilors %>% mutate(Sociocultural_threat = case_when(cultural_threat >= 1.410  ~ 'High',
                                                                     cultural_threat < 1.410  ~ 'Low'))
councilors$Sociocultural_threat <- factor(councilors$Sociocultural_threat, levels = c('High', 'Low'))
s1 <- plot(cj(councilors, f1, id = ~ResponseId, by = ~Sociocultural_threat, estimate = "mm"),
           group = "Sociocultural_threat", vline = 0.5)
s11 <- s1 + apatheme3 + ggtitle("Median point cut-off")

councilors <- councilors %>% mutate(Sociocultural_threat = case_when(cultural_threat >= 1.951  ~ 'High',
                                                                     cultural_threat < 1.951  ~ 'Low'))
councilors$Sociocultural_threat <- factor(councilors$Sociocultural_threat, levels = c('High', 'Low'))
s2 <- plot(cj(councilors, f1, id = ~ResponseId, by = ~Sociocultural_threat, estimate = "mm"),
           group = "Sociocultural_threat", vline = 0.5)
s22 <- s2 + apatheme3+ ggtitle("Mean point cut-off")

vTert = quantile(councilors$cultural_threat, c(0:3/3), na.rm = TRUE)
councilors$Sociocultural_threat = with(councilors, 
                           cut(cultural_threat, 
                               vTert, 
                               include.lowest = T, 
                               labels = c("Low", "Medium", "High")))
s3 <- plot(cj(councilors, f1, id = ~ResponseId, by = ~Sociocultural_threat, estimate = "mm"),
           group = "Sociocultural_threat", vline = 0.5)
s33 <- s3 + apatheme3+ ggtitle("Tertile cut-off")

grid.arrange(s11,s22,s33, ncol = 3, layout_matrix=rbind(c(1,1,2,2), c(NA, 3, 3, NA)))


# Figure D.11: Economic threat PCA sensitivity checks
summary(councilors$economic_threat)

councilors <- councilors %>% mutate(Economic_threat = case_when(economic_threat >= 3.525  ~ 'High',
                                                                     economic_threat < 3.525  ~ 'Low'))
councilors$Economic_threat <- factor(councilors$Economic_threat, levels = c('High', 'Low'))
s1 <- plot(cj(councilors, f1, id = ~ResponseId, by = ~Economic_threat, estimate = "mm"),
           group = "Economic_threat", vline = 0.5)
s11 <- s1 + apatheme3 + ggtitle("Median point cut-off")

councilors <- councilors %>% mutate(Economic_threat = case_when(economic_threat >= 3.778  ~ 'High',
                                                                     economic_threat < 3.778  ~ 'Low'))
councilors$Economic_threat <- factor(councilors$Economic_threat, levels = c('High', 'Low'))
s2 <- plot(cj(councilors, f1, id = ~ResponseId, by = ~Economic_threat, estimate = "mm"),
           group = "Economic_threat", vline = 0.5)
s22 <- s2 + apatheme3+ ggtitle("Mean point cut-off")

vTert = quantile(councilors$economic_threat, c(0:3/3), na.rm = TRUE)
councilors$Economic_threat = with(councilors, 
                                       cut(economic_threat, 
                                           vTert, 
                                           include.lowest = T, 
                                           labels = c("Low", "Medium", "High")))
s3 <- plot(cj(councilors, f1, id = ~ResponseId, by = ~Economic_threat, estimate = "mm"),
           group = "Economic_threat", vline = 0.5)
s33 <- s3 + apatheme3+ ggtitle("Tertile cut-off")

grid.arrange(s11,s22,s33, ncol = 3, layout_matrix=rbind(c(1,1,2,2), c(NA, 3, 3, NA)))


# Figure D.12: Ethnocentric values PCA sensitivity checks
summary(councilors$ethnocentric_pca)

councilors <- councilors %>% mutate(Ethnocentric_values = case_when(ethnocentric_pca >= 6.544  ~ 'High',
                                                                ethnocentric_pca < 6.544  ~ 'Low'))
councilors$Ethnocentric_values <- factor(councilors$Ethnocentric_values, levels = c('High', 'Low'))
s1 <- plot(cj(councilors, f1, id = ~ResponseId, by = ~Ethnocentric_values, estimate = "mm"),
           group = "Ethnocentric_values", vline = 0.5)
s11 <- s1 + apatheme3 + ggtitle("Median point cut-off")

councilors <- councilors %>% mutate(Ethnocentric_values = case_when(ethnocentric_pca >= 6.400  ~ 'High',
                                                                ethnocentric_pca < 6.400  ~ 'Low'))
councilors$Ethnocentric_values <- factor(councilors$Ethnocentric_values, levels = c('High', 'Low'))
s2 <- plot(cj(councilors, f1, id = ~ResponseId, by = ~Ethnocentric_values, estimate = "mm"),
           group = "Ethnocentric_values", vline = 0.5)
s22 <- s2 + apatheme3+ ggtitle("Mean point cut-off")

vTert = quantile(councilors$ethnocentric_pca, c(0:3/3), na.rm = TRUE)
councilors$Ethnocentric_values = with(councilors, 
                                  cut(ethnocentric_pca, 
                                      vTert, 
                                      include.lowest = T, 
                                      labels = c("Low", "Medium", "High")))
s3 <- plot(cj(councilors, f1, id = ~ResponseId, by = ~Ethnocentric_values, estimate = "mm"),
           group = "Ethnocentric_values", vline = 0.5)
s33 <- s3 + apatheme3+ ggtitle("Tertile cut-off")

grid.arrange(s11,s22,s33, ncol = 3, layout_matrix=rbind(c(1,1,2,2), c(NA, 3, 3, NA)))



# Figure D.13: Reported party affiliation vs. reported ideology of councilors

councilors2 <- councilors[!duplicated(councilors$ResponseId), ]

ReportedIdeology <- as.numeric(councilors2$Q16_1)


councilors2$Party[councilors2$kke==1]<-"KKE"
councilors2$Party[councilors2$mera25==1]<-"Mera25"
councilors2$Party[councilors2$pleusi==1]<-"Pleusi"
councilors2$Party[councilors2$syriza==1]<-"Syriza"
councilors2$Party[councilors2$kinal==1]<-"Kinal"
councilors2$Party[councilors2$nd==1]<-"New Democracy"
councilors2$Party[councilors2$greek_solution==1]<-"Party not indicated"
councilors2$Party[councilors2$golden_d==1]<-"Party not indicated"
councilors2$Party[is.na(councilors2$Party)]<-"Party not indicated"

councilors2$Party <- factor(councilors2$Party, levels = c("KKE", "Mera25", "Pleusi",
                                                          "Syriza", "Kinal", "New Democracy",
                                                          "Party not indicated"))

with( councilors2 , aggregate( Q16_1 , by=list(Party) , FUN=summary)  )

table(councilors2$Party)
table(councilors2$Party, councilors2$Q16_1)

# plot
p <- councilors2 %>%
  mutate(text = fct_reorder(Party, Q16_1)) %>%
  ggplot( aes(x=Q16_1, color=Party, fill=Party), na.rm = TRUE) +
  geom_histogram(alpha=0.6, binwidth = 1)+
  theme_ipsum() +
  theme(
    legend.position="none",
    panel.spacing = unit(0.2, "lines"),
    strip.text.x = element_text(size = 14)
  ) +
  xlab("Ideology") +
  facet_wrap(~text)
p


# Table D.14: Balance test (Covariates as dependent variables)

bal1 <- lm_robust(cultural_threat ~ Proximity + Public_goods + Size+ Type + Run_by, clusters = ResponseId,
                  data = councilors)

bal2 <- lm_robust(economic_threat ~ Proximity + Public_goods + Size+ Type + Run_by, clusters = ResponseId,
                  data = councilors)

bal3 <- lm_robust(ideology ~ Proximity + Public_goods + Size+ Type + Run_by, clusters = ResponseId,
                  data = councilors)

bal4 <- lm_robust(ethnocentric_pca ~ Proximity + Public_goods + Size+ Type + Run_by, clusters = ResponseId,
                  data = councilors)

bal5 <- lm_robust(treat1 ~ Proximity + Public_goods + Size+ Type + Run_by, clusters = ResponseId,
                  data = councilors)

texreg(list(bal1, bal2, bal3, bal4, bal5),
       custom.model.names = c("Sociocultural threat", "Economic threat",
                              "Ideology", "Ethnocentric values", "Active camp"), include.ci=FALSE)


